home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlae2.f < prev    next >
Encoding:
Text File  |  1997-06-25  |  3.4 KB  |  125 lines
  1.       SUBROUTINE DLAE2( A, B, C, RT1, RT2 )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       DOUBLE PRECISION   A, B, C, RT1, RT2
  10. *     ..
  11. *
  12. *  Purpose
  13. *  =======
  14. *
  15. *  DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
  16. *     [  A   B  ]
  17. *     [  B   C  ].
  18. *  On return, RT1 is the eigenvalue of larger absolute value, and RT2
  19. *  is the eigenvalue of smaller absolute value.
  20. *
  21. *  Arguments
  22. *  =========
  23. *
  24. *  A       (input) DOUBLE PRECISION
  25. *          The (1,1) element of the 2-by-2 matrix.
  26. *
  27. *  B       (input) DOUBLE PRECISION
  28. *          The (1,2) and (2,1) elements of the 2-by-2 matrix.
  29. *
  30. *  C       (input) DOUBLE PRECISION
  31. *          The (2,2) element of the 2-by-2 matrix.
  32. *
  33. *  RT1     (output) DOUBLE PRECISION
  34. *          The eigenvalue of larger absolute value.
  35. *
  36. *  RT2     (output) DOUBLE PRECISION
  37. *          The eigenvalue of smaller absolute value.
  38. *
  39. *  Further Details
  40. *  ===============
  41. *
  42. *  RT1 is accurate to a few ulps barring over/underflow.
  43. *
  44. *  RT2 may be inaccurate if there is massive cancellation in the
  45. *  determinant A*C-B*B; higher precision or correctly rounded or
  46. *  correctly truncated arithmetic would be needed to compute RT2
  47. *  accurately in all cases.
  48. *
  49. *  Overflow is possible only if RT1 is within a factor of 5 of overflow.
  50. *  Underflow is harmless if the input data is 0 or exceeds
  51. *     underflow_threshold / macheps.
  52. *
  53. * =====================================================================
  54. *
  55. *     .. Parameters ..
  56.       DOUBLE PRECISION   ONE
  57.       PARAMETER          ( ONE = 1.0D0 )
  58.       DOUBLE PRECISION   TWO
  59.       PARAMETER          ( TWO = 2.0D0 )
  60.       DOUBLE PRECISION   ZERO
  61.       PARAMETER          ( ZERO = 0.0D0 )
  62.       DOUBLE PRECISION   HALF
  63.       PARAMETER          ( HALF = 0.5D0 )
  64. *     ..
  65. *     .. Local Scalars ..
  66.       DOUBLE PRECISION   AB, ACMN, ACMX, ADF, DF, RT, SM, TB
  67. *     ..
  68. *     .. Intrinsic Functions ..
  69.       INTRINSIC          ABS, SQRT
  70. *     ..
  71. *     .. Executable Statements ..
  72. *
  73. *     Compute the eigenvalues
  74. *
  75.       SM = A + C
  76.       DF = A - C
  77.       ADF = ABS( DF )
  78.       TB = B + B
  79.       AB = ABS( TB )
  80.       IF( ABS( A ).GT.ABS( C ) ) THEN
  81.          ACMX = A
  82.          ACMN = C
  83.       ELSE
  84.          ACMX = C
  85.          ACMN = A
  86.       END IF
  87.       IF( ADF.GT.AB ) THEN
  88.          RT = ADF*SQRT( ONE+( AB / ADF )**2 )
  89.       ELSE IF( ADF.LT.AB ) THEN
  90.          RT = AB*SQRT( ONE+( ADF / AB )**2 )
  91.       ELSE
  92. *
  93. *        Includes case AB=ADF=0
  94. *
  95.          RT = AB*SQRT( TWO )
  96.       END IF
  97.       IF( SM.LT.ZERO ) THEN
  98.          RT1 = HALF*( SM-RT )
  99. *
  100. *        Order of execution important.
  101. *        To get fully accurate smaller eigenvalue,
  102. *        next line needs to be executed in higher precision.
  103. *
  104.          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  105.       ELSE IF( SM.GT.ZERO ) THEN
  106.          RT1 = HALF*( SM+RT )
  107. *
  108. *        Order of execution important.
  109. *        To get fully accurate smaller eigenvalue,
  110. *        next line needs to be executed in higher precision.
  111. *
  112.          RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
  113.       ELSE
  114. *
  115. *        Includes case RT1 = RT2 = 0
  116. *
  117.          RT1 = HALF*RT
  118.          RT2 = -HALF*RT
  119.       END IF
  120.       RETURN
  121. *
  122. *     End of DLAE2
  123. *
  124.       END
  125.